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Progress  has  been  made  on  the  following,  topics. 

1.  Averaged  fluid  equations:  the  LANS-o:  and  CLANS-a  equations 

2.  Simulations  of  the  compressible  CLANS-o  equations 

3.  Dynamic  a  formulation  and  computations  for  turbulent  flow 

4.  The  Helmholtz  deconvolution  model 

5.  Optimal  control  in  a  dynamic  environment 

6.  The  EPDiff  Equations 

7.  Flow  simulation  around  a  Micro  Air  Vehicle  (MAV) 

1  Averaged  Fluid  Equations:  the  LANS-ct  and  CLANS- 
a  Equations 

We  have  derived,  in  a  mathematically  precise  way,  Lagrangian  averaged  models  for  incom¬ 
pressible  and  compressible  flow.  This  work  has  now  appeared  in  the  SIAM  journal  Multiscale 
Modeling  and  Simulation ;  Bhat,  Fetecau,  Marsden,  Mohseni,  and  West  [2005].  That  work 
lays  important  foundations  for  our  computational  work  on  the  averaged  equations  as  a  com¬ 
putational  tool  for  shocks. 

The  LAE  -a  equations:  (due  to  Holm,  Marsden  and  Ratiu  in  1998)  are 

d 

— m  +  (u  ■  V)m  —  a2(Vu)T  •  A u  =  —  gradp, 

where  a  >  0  is  a  small  parameter  indicating  the  subgrid  scale,  m  =  (1  —  a2  A )u,  divrr  =  0, 
and  p  is  the  fluid  pressure. 

The  (much  simplified)  idea  of  the  derivation  of  EPDiff  and  LAE-a  is  as  follows: 

1.  average  Hamilton’s  principle  of  ideal  fluid  mechanics  over  trial  curves  forming  a  tube 
of  size  a. 

2.  Expand  in  terms  of  a  and  truncate  at  a2 

3.  Make  a  fluctuation  advection  hypothesis — a  Taylor  hypothesis  (the  fluctuations  are  Lie 
advected  by  the  main  flow) 
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4.  This  leads  to  the  action  function 


5.  (Compute  the  Euler-Poincare  equations  for  this  averaged  action)  =>  LAE-a  equations. 

While  this  basic  idea  for  the  derivation  has  been  around  for  some  time,  the  technical 
details  considerably  improves  earlier  ones,  even  in  the  incompressible  isotropic  case  and, 
moreover,  it  allowed  us  to  make  progress  on  the  compressible  case  as  well  (see  Bhat,  Fetecau, 
Marsden,  Mohseni,  and  West  [2005]).  A  simplified  version  of  what  we  are  doing  can  be  seen 
from  the  filtered  Lagrangian  for  Id  compressible  flow 

Kp>u)  =  J {\uv ~w(p)Sj  pdNx,  I1-1) 

where  v  =  (l  —  a2dxx)u. 

Obtain  the  equations  of  motion  from  this  filtered  Lagrangian  again  the  version  of  Euler- 
Poincare  theory  that  contains  provision  for  advected  quantities,  in  this  case  the  density. 

Sample  calculation  shown  in  Figure  1.1  on  equations  of  this  sort  show  that  a  controls  the 
width  of  the  shock,  independent  of  the  magnitude  of  (numerical)  hyperviscosity  that  is  used 
and  it  picks  up  the  correct  entropy  solution.  As  in  the  incompressible  turbulent  calculations, 
this  should  help  with  resolution  issues. 
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Figure  1.1:  Shock  regularization  by  Lagrangian  averaging  approach.  Left:  Burgers  shock 
regularization.  Right:  Normal  shock  regularization. 


2  Dynamic  LANS-ct  Model 

In  the  last  workshop  and  progress  report,  we  reported  on  progress  in  turbulent  flow  cal¬ 
culations  for  the  incompressible  case  in  which  averaged  isotropic  fluid  equations  accurately 
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compute  the  flow  at  scales  above  a  (see  Mohseni,  Kosovic,  Shkoller  and  Marsden  [2003]), 
but  does  not  have  to  resolve  scales  below  a.  “The  advantages  of  LES-type  models,  but  there 
are  well-defined  pde’s  that  keep  the  structure  of  fluid  mechanics.” 

We  developed  a  dynamic  procedure  for  the  Lagrangian  Averaged  Navier-Stokes-a  (LANS- 
a)  equations  where  the  variation  in  the  parameter  a  in  the  direction  of  anisotropy  is  deter¬ 
mined  in  a  self-consistent  way  from  data  contained  in  the  simulation  itself  (see  Zhao,  Mohseni 
and  Marsden  [2004];  Zhao  and  Mohseni  [2004];  Mohseni  and  Zhao  [2005]).  The  development 
of  the  dynamic  LANS-a  model  of  this  section  is  a  substantial  progress  in  application  of  the 
LANS-Q'  model  to  anisotropic  flows,  such  flow  around  an  airplane.  In  order  to  derive  this 
model,  the  incompressible  Navier-Stokes  equations  are  Helmholtz-filtered  at  the  grid  and  a 
test  filter  levels.  A  Germano  type  identity  is  derived  as  follows  by  comparing  the  filtered 
subgrid  scale  stress  terms  with  those  given  in  the  LANS-a  equations. 

Lij  Tij  'Tij  UilLj  V-ill j 

Where  (•)  stands  for  grid  filter  and  (•)  stands  for  test  filter. 

Assuming  constant  a  in  homogenous  directions  of  the  flow  and  averaging  in  these  direc¬ 
tions,  results  in  a  nonlinear  equation  for  the  parameter  a 

where  (5  is  the  aspect  ratio  of  test  and  grid  filter  width  and 

M  duj  duj  duk  duk  |  dv.j  duk 

13  dxk  dxk  dxt  dxj  dxk  dxj 

dui  duj  duk  duk  dui  duk 

yy.  = - ±  - - 1-  - - 

lJ  dxk  dxk  dxi  dxj  dxk  dxj 

Using  the  above  equations,  the  parameter  a  is  calculated  during  the  simulation  instead 
of  a  pre-defined  value.  Our  dynamic  model  is  initially  tested  in  forced  and  decaying  isotropic 
turbulent  flows  where  a  is  constant  in  space  but  it  is  allowed  to  vary  in  time.  It  is  observed 
that  by  using  the  dynamic  LANS-a  procedure  a  more  accurate  simulation  of  the  isotropic 
homogeneous  turbulence  is  achieved.  The  energy  spectra  and  the  total  kinetic  energy  decay 
are  captured  more  accurately  as  compared  with  the  LANS-a  simulations  using  a  fixed  a  (see 
Figure  2.1). 


a2  =  F{a) 


In  order  to  evaluate  the  applicability  of  the  dynamic  LANS-a  model  in  anisotropic  turbu¬ 
lence,  a  priori  test  of  a  turbulent  channel  flow  is  performed  (see  Figure  2.2).  It  is  found  that 
the  parameter  a  changes  in  the  wall  normal  direction.  Near  a  solid  wall,  the  length  scale  a 
is  seen  to  depend  on  the  distance  from  the  wall  with  a  vanishing  value  at  the  wall.  On  the 
other  hand,  away  from  the  wall,  where  the  turbulence  is  more  isotropic,  a  approaches  an  al¬ 
most  constant  value.  Furthermore,  the  behavior  of  the  subgrid  scale  stresses  in  the  near  wall 


3 


Figure  2.1:  Dynamic  LANS- a  simulation  for  isotropic  homogenous  turbulence.  Left:  Total 
kinetic  energy  decay.  Right:  Energy  spectra  snapshots. 

region  is  captured  accurately  by  the  dynamic  LANS- a  model.  The  dynamic  LANS-a  model 
has  the  potential  to  extend  the  applicability  of  the  LANS-a  equations  to  more  complicated 
anisotropic  flows. 


Figure  2.2:  A-priori  test  of  dynamic  LANS-a  model  for  turbulent  channel  flow.  Left:  Vari¬ 
ations  of  a  in  wall  normal  direction.  Right:  Sub- grid  stress  distribution  in  wall  normal 
direction. 


3  Simulations  of  the  Compressible  LANS-a  Equations 

Numerical  Solution  and  Validation  of  CLANS-a:  In  order  to  test  and  validate  the 
performance  of  the  CLANS-a  equations,  one  first  needs  to  prepare  a  database  for  direct 
numerical  simulation  of  compressible  isotropic  turbulence  and  shock  turbulence  interaction 
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using  the  Navier-Stokes  equations.  A  parallel  three-dimensional  code,  capable  of  simulating 
compressible  isotropic  turbulence  and  shock-turbulence  interaction  was  developed  using  a 
combination  of  WENO  (Weighted  Essentially  Non-Oscillating)  and  compact  scheme.  This 
scheme  has  fifth  order  accuracy  in  space  and  it  employs  a  third  order  Runge-Kutta  scheme 
for  time  integration.  WENO  scheme  is  an  efficient  technique  in  shock  capturing  problems 
and  is  particularly  effective  in  simulating  shock-turbulence  interaction.  We  have  tested 
and  validated  our  DNS  solver  by  comparing  the  results  with  the  ones  obtained  by  Pullin, 
Samtaney  and  Kosovic  [2001].  Figures  3.1  show  the  comparisons  of  the  total  kinetic  energy 
decay  and  time  evolution  of  the  energy  spectra  in  a  decaying  turbulence  case. 


Figure  3.1:  Direct  numerical  simulation  of  compressible  turbulence.  Left:  Total  kinetic 
energy  decay.  Right:  Energy  spectra  snapshots  and  comparision  with  k  ~  —  5/3  law. 


1-D  Inviscid  Shock  Regularization  and  Shock- Turbulence  Interaction  Simulation 
Using  the  Averaged  Euler  Equations:  A  preliminary  study  of  inviscid  shock  regular¬ 
ization  by  averaging  the  convective  velocity  in  Euler  equations  was  performed  on  a  1-D  shock 
tube  problem.  It  is  observed  that  the  shock  in  pressure  are  nicely  regularized  (see  Figure  3.2) 
by  replacing  the  convective  term  with  a  Helmholtz  filtered  velocity  u  =  (1  —  a2A)_1u.  In  this 
inviscid  shock  regularization,  the  shock  thickness  is  controlled  by  averaging  length  scale  a. 
We  also  studied  one- dimensional  shock-turbulence  interaction  by  this  averaged  Euler  equa¬ 
tions.  It  is  found  that  pressure  fluctuation  after  the  shock  can  be  smoothed  by  increasing 
the  a  value  in  the  averaged  Euler  equation,  (see  Figure  3.2) 


Large  Eddy  Simulation  of  Compressible  Navier-Stokes  Equations:  Based  on  our 
DNS  code,  we  also  developed  a  code  for  large  eddy  simulation  of  compressible  turbulence  to 
be  used  for  testing  the  accuracy  and  computational  cost  of  our  3-D  CLANS-a  model.  This 
code  will  be  used  to  evaluate  the  performance  of  our  CLANS-a  solver.  We  have  implemented 
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Figure  3.2:  Shock  regularization  and  shock-turbulence  interaction.  Left:  Pressure  distribu¬ 
tion  in  a  shock  tube.  Right:  One-dimensional  shock-turbulence  interaction. 

a  Smagorinsky  and  dynamic  Smagorinsky  model  for  sub-grid  scale  modeling.  Figure  3.3  and 
Figure  3.4  show  the  LES  and  dynamic  LES  simulation  results  as  well  as  the  comparisons 
with  DNS  results  for  isotropic  homogenous  turbulence. 


Figure  3.3:  Large-eddy  simulation  of  compressible  turbulence.  Left:  Total  kinetic  energy 
decay.  Right:  Energy  spectra  snapshots. 


4  Helmholtz  Deconvolution  Model 


We  have  already  shown  that  the  LANS-a  model  can  accurately  capture  turbulence  scales 
larger  than  the  scale  a  Mohseni,  Kosovic,  Shkoller  and  Marsden  [2003].  Since  the  Helmholtz 
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Figure  3.4:  Large-eddy  simulation  of  compressible  turbulence.  Left:  Total  kinetic  energy 
decay.  Right:  Energy  spectra  snapshots. 


filter  is  inver table,  the  turbulence  scales  below  a  may  be  recovered  by  a  deconvolution 
operation.  This  leads  us  to  derive  a  deconvolution  model  with  Helmholtz  filter. 

This  idea  was  initiated  by  A.  Jameson  (Stanford  University)  based  on  our  Lagrangian 
averaging  technique.  A  collaboration  with  Jameson  started  in  the  last  AFOSR  PI  meeting 
in  Dayton  Ohio.  Following  the  meeting,  we  developed  a  Helmholtz  deconvolution  model  (see 
Zhao,  Mohseni  and  Jameson  [2004]).  By  substituting  Helmholtz  decomposed  terms  U{  = 
(1  -  a2  A )Hi  into  the  incompressible  Navier-Stokes  equations  and  rearranging  the  equations, 

we  can  write  _  _  _ 

dui  duiUj  _  dp  d2u,i  dr^ 

^  dt  ^  dxj  dxi  +  ^  dxjdxj  ^  dxj 


where 


Tij  =  a2(l  -  a2 A)  1 


„  dUi  diij  9  d2Ui  d2Uj 

2 — - — 1  +  a2 - - - J— 

dxk  dxk  dxkdxk  dxtdxt 


In  order  to  avoid  inverting  the  Helmholtz  operator,  one  can  expand  the  operator  as 


(1  -  a2A)  1  =  1  +  a2 A  +  a4 A2  -| - 


By  retaining  terms  up  to  the  fourth  power  of  a,  the  approximate  virtual  stress  tensor  assumes 
the  following  form 


=  2a2^^  +  a4 
dxk  axk 


+  AuiAuj 


We  performed  some  numerical  simulations  with  this  model  for  isotropic  homogenous 
turbulence.  Figure  4.1  shows  our  preliminary  results  with  spectral  vanishing  viscosity  (SVV). 
Our  Helmholtz  deconvolution  model  provides  an  accurate  sub-grid  scale  stress  model  for  the 
decaying  isotropic  homogenous  turbulence  simulation  in  this  study.  We  are  in  the  process  of 
applying  this  model  to  more  sever  test  case. 
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Figure  4.1:  Helmholtz  deconvolution  model  in  Zhao,  Mohseni  and  Jameson  [2004]  for 
isotropic  homogenous  turbulence  simulation.  Left:  Total  kinetic  energy  decay.  Right:  En¬ 
ergy  spectra  snapshots. 

5  Optimal  Control  in  a  Dynamic  Environment 

The  objective  of  this  project  is  to  do  an  optimization  of  underwater  gliders  in  a  time  depen¬ 
dent  ocean  environment.  This  study  was  motivated  by  the  AOSN-II  project  (Autonomous 
Ocean  Sampling  Network).  The  results  are  given  in  Inane,  Shadden,  and  Marsden  [2005]. 

The  problem  involves  the  control  of  a  simplified  model  of  an  underwater  glider  in  Mon¬ 
terey  Bay  moving  in  ocean  currents  whose  velocities  can  be  comparable  to  that  of  the  glider 
itself. 


The  objective  is  to  get  from  point  A  (the  top  “cross”  in  the  figure  below)  to  point  B  (the 
bottom  “cross”),  satisfying  model  equations  of  motion,  and  at  the  same  time,  optimizing  a 
cost  function.  In  this  case,  the  cost  function  was  a  weighted  sum  of  a  temporal  cost  and 
an  energy  cost.  We  then  compare  the  resulting  path  with  Lagrangian  coherent  structures 
(LCS),  which  are  regions  of  expanding  (or  attracting)  material  lines  and  play  the  role,  in  the 
time  dependent  context,  of  stable  and  unstable  manifolds.  We  wanted  to  test  a  conjecture, 
motivated  by  the  use  of  natural  dynamics  in  other  areas,  that  optimal  trajectories  are  in 
fact,  linked  to  the  use  of  LCS  pathways  in  the  ocean.  In  this  study,  this  conjecture  was 
verified. 
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The  technique  involves  several  key  steps.  First  of  all,  we  gather  data  using  the  Harvard 
HOPS  model  for  data  assimilated  from  Monterey  Bay  last  summer  in  August,  2003.  We 
use  this  data  and  the  software  MANGEN  (Manifold  Generation — Lekien  and  Couliette)  to 
compute  Lagrangian  coherent  structures  for  the  velocity  field  so  collected.  With  the  same 
velocity  field,  we  also  use  the  software  NTG  (Nonlinear  Trajectory  Generation)  to  compute 
the  optimal  trajectory.  Finally,  we  compare  these  two  independent  calculations  and  see  if 
there  is  agreement.  Indeed,  the  agreement  is  remarkably  good,  indicating  that  one  can  use 
the  LCS  pathways  for  optimal  planning  purposes. 


36  » 
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Longitude  (degrees) 


Shane  Ross  also  helped  us  transition  some  of  the  Mangen  capabilities  from  the  fluid 
context  to  that  of  transport  and  mixing  in  solar  system  dynamics. 


6  The  EPDiff  Equations 


The  work  described  in  this  section  is  given  in  Holm  and  Marsden  [2004].  The  n-dimensional 
EPDiff  equations  are  as  follows 


U  .  .  _  rv 

— m  +  •  Vm,  +  .  Vu  ^  •  mj 

convection  stretching 


+  m  (div  u ) 

' - V - ' 

expansion 


0,  the 


where  m  =  (1  —  ot2A)u. 

These  equations  fit  into  the  general  theme  of  averaged  equations  in  fluid  mechanics  (Euler- 
Poincare  equations;  eg,  variational,  Hamiltonian  structure,  Kelvin  theorem,  etc.).  In  fact, 
for  n  =  2,  they  serve  as  a  model  for  2-dimensional  waves 

These  equations  have  interesting  special  interface  solutions — one  of  our  main  contribu¬ 
tions  is  to  reveal  the  rich  geometric  structure  about  singular  solutions  by  showing  that  the 
solution  Ansatz  is  given  by  a  momentum  map  relation. 
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The  computation  of  solutions  were  done  by  Martin  Staley  at  Los  Alamos — and  makes 
use  of  mimetic  differencing  methods;  this  approach  is  important  to  get  good  resolution  of 
the  structures  and  is  related  to  our  idea  to  eventually  implement  variational  integrators  in 
fluid  problems. 

The  equations  in  the  1-d  case  are  the  Camassa-Holm  equations  for  shallow  water  waves 
and  the  solutions  are,  as  is  well-known,  dominated  by  peakons,  as  in  the  figure. 


The  2-d  case,  on  the  other  hand,  is  related  to  two  dimensional  waves.  The  figure  shows 
interface  solutions  of  EPDiff  compared  to  synthetic  aperture  radar  observations  (from  the 
Space  Shuttle)  of  internal  waves  in  the  South  China  Sea.  It  is  clear  that  these  special  interface 
solutions  of  the  EPDiff  equation  are  important. 
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In  summary,  what  we  did  was: 

1.  Put  the  problem  into  the  Euler-Poincare  framework  (geometric  fluid  mechanics) 

2.  Discover  the  geometry  of  the  interface  solutions  (eg,  that  they  also  evolve  by  a  me¬ 
chanical  structure,  have  their  own  conservation  laws,  etc). 


7  Flow  Simulation  around  a  Micro  Air  Vehicle  (MAY) 


A  long  term  goal  of  our  group  is  to  implement  the  LANS-a  model  in  simulating  turbulent 
flows  in  complex  geometry.  The  dynamic  LANS-a  model  reported  in  section  2  is  moving  us 
in  that  direction.  The  dynamic  LANS-a  model  allows  us  to  extend  the  LANS-a  calculations 
to  anisotropic  flows  e.g.,  close  to  a  wall.  We  have  started  to  develope  our  computational 
capability  for  simulating  flows  in  complex  geometry.  This  section  summarize  this  effort  (also 
see  Mohseni  et  al  [2004]). 

A  micro  air  vehicle  (MAV)  with  fixed  flexible  wings  was  recently  designed,  built,  and 
tested  in  Mohseni ’s  group  (see  Figure  7.1).  The  MAV  is  designed  to  serve  as  a  platform 
for  characterizing  the  flight  performance  of  different  low  Reynolds  number  wing  designs  by 
measuring  the  in-flight  accelerations  due  to  control  inputs  and  gusts.  The  platform  houses 
an  onboard  flight  computer  and  communications  system  that  can  process  sensor  information, 
transmit  data  to  a  ground  station,  and  control  the  aircraft  via  pilot  input  from  the  ground 
station.  With  light  weight  and  small  wing  aspect  ratio,  the  Colorado  MAV  is  susceptible  to 
atmospheric  turbulence  and  gust  wind  conditions.  For  optimal  design  of  such  a  vehicle  it  is 
necessary  to  calculate  the  aerodynamic  characteristice  of  the  MAV. 


Figure  7.1:  Colorado  MAV,  a  flexible  wing  micro  air  vehicle  designed,  fabricated,  and  tested 
at  the  University  of  Colorado,  (a)  CAD  model  assembly,  (b)  CAD  model,  (c)  CMAV  proto¬ 
type. 

A  domain  decomposition  based  parallel  three-dimensional  compressible  Euler  and  Navier- 
Stokes  CFD  solver  is  developed  which  combines  finite  volume  and  finite  element  discretiza¬ 
tions  on  unstructured  tetrahedral  meshes.  It  features  an  upwind  scheme  using  a  piecewise 
linear  reconstruction  of  the  flow  variables  in  each  control  volume  for  the  convective  term, 
and  a  PI  finite  element  Galerkin  approximation  for  the  diffusive  term. 


11 


Starting  from  the  CAD  model  shown  in  Figure  7.1(b),  an  anisotropic  grid  was  generated 
with  459,000  points  and  2.1  million  tetrahedra.  About  10  layers  of  stretched  elements  are 
located  close  to  the  MAV  in  order  to  capture  accurately  the  boundary  layers.  The  average 
thickness  of  these  stretched  elements  is  equal  to  10  mm.  Figure  7.2  displays  several  views  of 
the  surface  mesh  as  well  as  a  cut  along  the  symmetry  plane  of  the  MAV.  The  refinement  on 
the  leading  and  trailing  edges  of  the  wing  as  well  as  on  the  engine  can  be  observed  in  these 
figures. 


Figure  7.2:  Computational  grid  around  the  CMAV  from  Mohseni  et  al  [2004]. 


The  flow  is  assumed  to  be  viscous  and  the  freestream  velocity  and  the  angle  of  attack  are 
first  set  respectively  to  14  m/s  and  15  deg.  Figure  7.3  shows  the  streamlines  obtained  for 
these  conditions.  The  freestream  velocity  and  the  angle  of  attack  are  then  varied.  Figure 
7.4  shows  the  aerodynamic  efficiency  curve,  where  a  maximum  glide  ratio  of  3.8  is  observed. 
Similar  glide  ratio  was  measured  by  Waszak  and  Jenkins  [2001]  for  a  flexible  wing  MAV. 


Figure  7.3:  Streamlines  around  the  CMAV  at  15°  AoA.  From  Mohseni  et  al  [2004] 


Increase  in  flight  Reynolds  number  would  increase  the  lift  to  drag  ratio. 

Since  the  dynamic  LANS-a  model  of  section  2  preserve  the  geometrical  structure  of  the 
governing  equations  (e.g,  Kelvin  circulation)  we  expect  that  this  model  will  perform  much 
better  in  turbulence  simulation  around  MAVs.  We  are  now  in  the  process  of  implementing  a 
dynamic  LANS-a  model  for  the  turbulence  simulation  around  the  Colorado  MAV.  This  step 
will  move  us  closer  to  practical  application  of  the  LANS-a  modeling  in  problems  of  interest 
in  Aerospace  industry. 
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